package fem2.material;

import inf.jlinalg.SolveFailedException;
import inf.text.ArrayFormat;

import org.junit.Test;

import fem2.Constraint;
import fem2.Element;
import fem2.Face;
import fem2.MaterialModel;
import fem2.Mesh;
import fem2.MeshUtilities;
import fem2.Model;
import fem2.Node;
import fem2.PatchTestQ4;
import fem2.PatchTestT3;
import fem2.analysis.Analysis;
import fem2.analysis.NonlinearStaticManualAnalysis;
import fem2.element.BBarStructuralElement;
import fem2.element.StructuralElement;
import fem2.enu.EchoLevelType;
import fem2.enu.IntegrationType;
import fem2.enu.State;
import fem2.material.DruckerPragerMM;
import fem2.material.FluencyCriteria;
import fem2.material.LinearIsotropicHardeningRule;

//@RunWith(Parameterized.class)
public class PatchTestDruckerPrager {

	private Model model;

	// @Parameters
	// public static Collection<Object[]> parameters() {
	// List<Object[]> params = new LinkedList<Object[]>();
	//
	// params.add(new Object[] { createModel(createMeshQ4(), 0) });
	// params.add(new Object[] { createModel(createMeshT3(), 0) });
	// params.add(new Object[] { createModel(createMeshQ4(), 1) });
	// params.add(new Object[] { createModel(createMeshT3(), 1) });
	// // params.add(new Object[] { createModel(createMeshQ4(), 2) }); // fail
	// // params.add(new Object[] { createModel(createMeshT3(), 2) }); // fail
	// return params;
	// }

	// public PatchTestDruckerPrager(Model m) {
	// // model = m;
	// model = createModel(createMeshT3(), 0);
	// }

	public static MaterialModel createMaterial() {
		double E = 1e7;
		double nu = 0.48;
		double thickness = 1;
		double rho = 0;
		State ss = State.PLANE_STRAIN;
		// return new StVenantKirchhoffMM(E, nu, thickness, rho, ss);

		double phi = (20.0 / 180) * Math.PI;
		double psi = phi;
		/*
		 * inner cone
		 */
		double eta = 6 * Math.sin(phi) / (Math.sqrt(3) * (3 + Math.sin(phi)));
		double xi = 6 * Math.cos(phi) / (Math.sqrt(3) * (3 + Math.sin(phi)));
		double etabar = 6 * Math.sin(psi) / (Math.sqrt(3) * (3 + Math.sin(psi)));

		double c = 490;
		double H = 0;

		FluencyCriteria hardeningLaw = new LinearIsotropicHardeningRule(c, H);

		return new DruckerPragerMM(E, nu, thickness, rho, ss, eta, etabar, xi, hardeningLaw);
		// return new DruckerPragerMM3(E, nu, thickness, rho, ss, eta, etabar,
		// xi, hardeningLaw);
	}

	public static Model createModel(Mesh mesh, int elementType) {
		Model M = new Model(mesh);
		MaterialModel mm = createMaterial();

		for (int i = 0; i < mesh.countFaces(); i++) {
			Element e;

			switch (elementType) {
			case 0:
				e = new StructuralElement(mesh.getFace(i), mm);
				((StructuralElement) e).setIntegrationType(IntegrationType.REDUCED_INTEGRATION);
				break;
			case 1:
				e = new BBarStructuralElement(mesh.getFace(i), mm);
				break;
			// case 2:
			// e = new EASStructuralElement(mesh.getFace(i), mm,
			// EASModeType.Q1E4);
			// break;
			default:
				throw new Error("invalid element type");
			}

			M.addElement(e);

			e.setEchoLevel(EchoLevelType.DEBUG);
		}

		return M;
	}

	@Test
	public void testSolve1() throws SolveFailedException {
		model = createModel(createMeshQ8(), 0);
		Mesh mesh = model.getMesh();

		double du = -0.5;
		Constraint cu = new Constraint(true, false);
		cu.setValue(1, du);

		mesh.addConstraint(cu, MeshUtilities.seekNodesOnSurface(mesh, 0.0, 1.0, -1.0));

		double[] lambda = new double[] { 0.5, 1.0 };
		Analysis an = new NonlinearStaticManualAnalysis(model, lambda);
		an.getStrategy().setEchoLevel(EchoLevelType.OUTPUT);
		an.run();

		for (int i = 0; i < 4; i++) {
			double[] uHatLocal = model.getMesh().getNode(i).getUHatLocal();
			System.out.println(ArrayFormat.format(uHatLocal));
		}

		// double[][] Uh = new double[][] { { 0, 0 }, { 2.0, 0 }, { 2.0,
		// -0.967228028754730 },
		// { 0, -0.967228028754730 } };
		// for (int i = 0; i < 4; i++) {
		// double[] uHatLocal = model.getMesh().getNode(i).getUHatLocal();
		// Assert.assertEquals(Uh[i][0], uHatLocal[0], 1e-10);
		// Assert.assertEquals(Uh[i][1], uHatLocal[1], 1e-10);
		// }

		// IVector epsilon = ((StructuralElement)
		// m.getElement(0)).computeTotalStrain(0);
		// IVector sigma = ((StructuralElement) m.getElement(0)).getStress(0);
		// System.out.println("strain = " + epsilon);
		// System.out.println("stress = " + sigma);
	}

	public static Mesh createMeshQ4() {
		return new PatchTestQ4(0).createMesh();
	}

	public static Mesh createMeshT3() {
		return new PatchTestT3(0).createMesh();
	}

	public static Mesh createMeshQ8() {
		Mesh mesh = new Mesh();

		double l = 1;
		double h = 1;
		Node n1 = mesh.addNode(new Node(0, 0));
		Node n2 = mesh.addNode(new Node(l, 0));
		Node n3 = mesh.addNode(new Node(l, h));
		Node n4 = mesh.addNode(new Node(0, h));
		Node n5 = mesh.addNode(new Node(l / 2, 0));
		Node n6 = mesh.addNode(new Node(l, h / 2));
		Node n7 = mesh.addNode(new Node(l / 2, h));
		Node n8 = mesh.addNode(new Node(0, h / 2));

		mesh.addFace(new Face(n1, n2, n3, n4, n5, n6, n7, n8));

		mesh.setNumberOfFields(2);

		Constraint cx = new Constraint(false, true);
		Constraint cy = new Constraint(true, false);

		n1.addConstraint(cx);
		n1.addConstraint(cy);
		n5.addConstraint(cy);
		n2.addConstraint(cy);
		n8.addConstraint(cx);
		n4.addConstraint(cx);

		return mesh;
	}
}
